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Abstract 


A  numerical  method  for  the  computation  of  the  singular  behavior  of  the 
solution  of  the  Laplace  equation  is  proposed.  It  is  shown  that  the  accuracy 
of  the  computed  stress  intensity  factor  by  the  h,p  and  h-p  version  of  the 
finite  element  method  has  the  same  order  as  the  squau'e  of  the  error  of  the 
solution  measured  in  the  energy  norm.  Numerical  examples  au'e  given. 


1. 


Introduction 


It  is  well  known  that  the  solution  of  elliptic  pairtial  differential 
equations  is  singulzu'  in  the  neighborhood  of  the  edges  and  the  vertices  of  the 

3 

domain  of  definition  llclR  .  The  chau'acter  of  the  solution  can  be  described 
by  the  decomposition  of  the  solution  in  singuleu'  smd  regular  parts  (see  e.g. 
[D],[G1],  [G2],  [KO],  [P],  [PSl],  (PS2]).  Singulau'  behavior  of  the  solution 
is  of  large  importsince  in  many  applications.  It  is,  for  example,  directly 
related  to  the  problems  of  fracture  mechanics.  Hence  the  numerical  deter¬ 
mination  of  the  parameters  of  the  singulair  behavior  of  the  solution  is  of 
great  interest,  for  exaimple,  in  problems  of  structural  mechanics. 

The  major  tool  of  computational  structure  mechanics  is  the  finite  element 
method.  In  the  3  dimensional  analysis,  one  of  the  most  laborious  parts  of  the 
finite  element  computation  is  the  mesh  generation.  Hence  the  method  for  the 
determination  of  the  singular  parts  of  the  solution  should  be  fully  integrated 
with  the  data  and  algorithm  used  for  solving  the  boundary  value  problem  of  the 
partial  differential  generation  of  interest.  In  the  engineering  and 
mathematical  literature,  meuiy  methods  for  the  computation  of  stress  intensity 
factors  were  proposed.  Most  methods  use  different  approaches  for  the  approxi¬ 
mation  of  the  solution,  the  approximation  of  singularity  functions  (and 
adjoint  singularity  functions),  and  the  extraction  of  the  intensity  factors. 

Our  analysis  addresses  the  error  of  the  complete  approach,  i.e.,  it 
includes  the  error  of  the  finite  element  approximation,  the  error  of  the 
computed  singularity  function,  and  the  extraction  of  the  vertex  intensity 
factors.  The  method  analyzed  here  is  partially  related  to  the  ideas  in  [LN]. 
There,  however,  a  different  eigenvalue  problem  was  used  eind  no  error  analysis 
was  performed. 
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In  this  paper  we  propose  and  ainalyze  such  method  for  the  cheu'acteris.tiza- 
tlon  of  the  singulairity  in  the  neighborhood  of  the  vertex  of  the  domain 
ncR.  We  restrict  ourselves  here  to  the  Laplace  equation  and  polyhedral 
domains  only.  This  paper  is  the  first  in  a  series.  The  other  papers  will 
deal  with  elasticity  problems,  which  are  of  especially  l2u*ge  interest  in 
engineering.  This  method  was  implemented  in  the  prograim  STRIPE  [S]  and  a 
survey  of  the  results  in  connection  with  the  einalysis  of  complex  airplaine 
structures  is  given  in  [A]. 

In  the  neighborhood  of  a  vertex  the  solution  u  of  the  boundary  value 
problem  can  be  written  in  the  form  (see  (2.11)) 

(1.1)  u(x)  =  u_(x)  +  EC.S.(x) 

u  J  J 


where  C.  deptends  (globally)  on  the  solution  auid  S.(x)  depends  on  the 
d  d 

geometry  only.  The  so-called  stress  intensity  functions  Sj(x)  as  well  as  the 
so-called  stress  intensity  factors  can  be  computed  only  numerically.  The 

function  u^  in  (1.1)  veinishes  faster  towards  the  vertex  thaui  the  functions 
Sj(x)  (see  e.g.  (2.11b)  for  exact  formulation). 

The  solution  u(x)  is  computed  approximately  by  the  finite  element 


method.  The  error  of  the  finite  element  solution  u^  satisfies  aji  asymptotic 
convergence  estimate  of  the  form 


(1.2)  ||u  -u||  <cF(N(q)) 

where  ||•||  is  typically  the  energy  norm,  N(q)  is  the  number  of  used  degrees 

of  freedom  and  F((;)  is  a  decreasing  function  depending  on  the  used  method, 

e.g.,  h,  p  or  h-p  version.  In  practice  we  usually  see  in  (1.2)  approximate 

equality,  i.e.  «  instead  <  .  We  will  show  in  this  paper  that  the  stress 

2 

intensity  factors  cam  be  computed  with  the  accuracy  F(N(q))  ,  i.e.,  denoting 
by  the  finite  element  approximation  of  we  get 
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-Cjl  <C  F(N(q))^ 

where  C  is  independent  of  q. 

This  is  one  of  the  major  results  of  this  paper. 

Section  2  gives  the  formulation  of  the  problem.  Section  3  introduces  a 
Steklov  problem  and  shows  that  the  functions  Sj  in  (1.1)  au'e  solutions  of  the 
problem.  In  Section  4  we  derive  a  formula  for  the  extraction  of  the  exact 
stress  intensity  factors  from  the  exact  solution  u. 

Section  5  elaborates  on  the  finite  element  method  and  formulates  some 
assumption  about  the  meshes  used.  It  shows  that  these  assumptions  are  valid 
for  the  standard  h,  p  aind  h-p  versions  of  the  finite  element  method. 

Section  6  elaborates  on  the  numerical  computation  of  the  function  S^ 
in  (1.1)  auid  gives  the  estimates  of  the  error. 

Section  7  describes  the  numerical  computation  of  the  stress  intensity 
factors  Cj  and  proves  the  error  estimate.  It  shows  that  the  accuracy  oi’ 

is  of  the  same  order  as  the  square  of  the  error  of  the  finite  element 
solution  u^  when  measured  in  the  energy  norm. 

Section  8  presents  an  illustrative  example  computed  by  the  p-version  of 
the  finite  element  method  implemented  in  the  progreun  STRIPE. 
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2.  Formulation  of  the  problem 

I 

Let  ncR  be  a  polyhedron  with  the  boundary  an  =  U  F,  =  F,  where  F 

^=1  ^  ^ 

are  the  planar  (open)  faces  of  an.  By  Vj,,i  =  we  denote  the 

vertices  of  n  and  will  assume  that  the  vertex  V  =  is  located  in  the 
origin  (i.e.  =  0).  By  =  1,2, •••,n  we  denote  the  (open)  edges  of  n. 

E^,i  =  <,  •••,nQ,  are  all  the  edges  containing  the  vertex  V^.  Let 
5  =  min{dist(Ej^,  Vq)}  ,  £=  nQ  +  l,*-*,n.  Further  by  Q^,p>0  we  denote  the 
open  ball  with  the  center  in  and  radius  p.  Set  R  =  min(l,^  6)  then 
obviously  ~  ^R^^  where  K  is  the  infinite  cone  coinciding  with 

n  in  the  neighborhood  of  V..  For  any  p<5  we  denote  n  =  Q  nn,  = 

0  p  p  p 

aQnn,  F^  =  FnQ^.  Further  let  F  =  F^  u  F^^  where  =  U  F^^  =  U  F^, 

DaN  =  0  be  the  Dirichlet  and  Neumann  part  of  F  respectively  auid 

F  _  =  F^nQ  ,  F  =  F.,  aQ  . 
p,D  D  p  p,N  N  p 

We  will  be  interested  in  the  (weak)  solution  of  the  boundary  value 
problem 


(2.  la) 


-Au  =  g 


in  n 


(2.1b)  u  =  0  on  Fjj 

'  =  •  Is  '  8n  ''n' 

Denote  by  H^(n)  the  usual  Sobolev  space,  and  let 
(2.2a)  H^(n)  =  (u€  H^(Jl)  |u  =  0  on  F^^}  if  F^  * « 

(2.2b)  H^(n)  =  {u€H^(n)|  J  u  dx  =  0>  if  F^^  =  ^. 

n 

Further  let 

(2.3)  Bjj(u,v)  =  J  7u  •  7v  dx 

n 
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be  the  bilinear  form  defined  on  x  H^{n)  and 

(2.4)  II^IIeCQ)  =  (Bjj(u.u))'^^ 

Obviously  ||  •  is  equivalent  with  the  standau'd  Sobolev  norm 

(if  rp*0). 


II  •  II  1 
H^(n) 


Later  we  will  also  use  the  notation 

idin  )  =  {u€ H^(n  ) lu  =  0  on  r  if  r 

up  p  P.D  P.D  ^ 

Hj,(n  )  =  {u€H^(n  Ilf  u  dx  =  0}  if  r  ^  = 

up  P  J  P.D 

n 

p 

The  meaning  of  (u.v)  amd  ll^ll£(jj  \  is  obvious. 

Now  the  weak  formulation  of  the  problem  (2.1)-(2.3)  is:  Find  u6H^(n) 
such  that 


Remark  2.1.  We  assumed  that  u  =  0  on  for  simplicity  only.  The 

aissumption  (2.6)  (together  with  u  =  0  on  r„  )  is  more  essential.  We 

U 

will  briefly  comment  on  it  in  Section  7.  D 
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Together  with  the  usual  c^tesian  coordinates  x  =  (x^.x^.x^),  we  also 
will  use  in  the  neighborhood  of  (specifically  in  0^)  the  spherical 
(polau*)  coordinate  system,  (r.Q.y)  centered  in  the  vertex  V^, 


Let  y 

=  {S> 

be  the  set  of  all 

functions 

SeH^I 

(2.7) 

S  =  r^w(0,  v>) ,  A  >  - 

1 

2 

ajid  satisfying  (in 

the  weaJc  form)  the 

equations 

(2.8a) 

-AS  =  0 

in 

«R 

(2.8b) 

S  =  0 

on 

^R,D 

(2.8c) 

on 

^R,N‘ 

In  the  case 

when  all  the  faces  r.,J 

=  1,  •  • • , n 

which 

the  Neumann  part  of  the  boundary  (i.e.  r„  ^  =  ^),  the  constant  function  also 

K,  U 

satisfies  (2.8)  but  we  do  not  include  it  in  !f. 

Below  we  will  see  that  the  set  ^  is  not  empty  ajid  is  denumerable,  i.e. 
y  =  with  Sj  =  r^w^(0,y),  Aj  e  R.  We  will  assume  that  the  Sj 

are  ordered  such  that  The  functions  Se/  will  be  called 

singularity  functions. 

Since 

3S  .  -Ic 

the  function  S  satisfies  the  equation 


(2.9) 

where 


(2. 10) 


bj^(u,  v) 


=  R 


-1 


uv  ds 


R 
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It  is  well  known  (see  e.g.  (DJ  (G1 ] , [G2] , [KO] , [P] )  that  the  solution  u 
of  (2.1)  (or  equivalently  (2.5)>  admits  (under  aLssumption  (2.6)  for  any  s>0 
the  following  decomposition  on  £lp^: 

(2.na)  u  =  uj,.  Y. 


wnere 


(2. 11b) 


Aj+l/2<s 


1  7Uq  I  ^r  ^®dv  <  00. 


Here  (2.11b)  relates  to  the  seminorm  of  the  weighed  Sobolev  space  H^^(nj^). 
(2.11b)  implies  that  u^  is  in  modulo  constants.  If  U* 

then  UQ€H^^(Qp). 

We  will  obtain  (2.11)  as  a  consequence  of  the  results  in  Section  3  aind  4. 

Equation  (2.11)  shows  that  the  behavior  of  the  solution  near  the  vertex  is 

determined  by  the  singularity  functions  S.. 

\) 

Let 

w  =  >  C.S., 

s  J  J 

A .+-Ss 

J  2 

then  the  relative  error  in  the  energy  norm  between  u  auid  w^  goes  to  zero 

for  r — >0,  because  of  (2.11b).  Therefore  w^  is  a  good  approximation  of  u 

in  H^(£2j^)  if  r  is  sufficiently  small. 

The  numbers  C  in  (2.11)  are  called  the  vertex  intensity  factors. 

With  each  singularity  function  S.  of  the  form  (2.7)  we  will  associate 

0 

the  adjoint  singularity  function 


(2. 12) 


S_  =  r~^  '^Jw.O.ip) 
J  J 


It  is  easy  to  check  that  S  .  satisfies  (2.8  abc),  but  S 

~J  -j  D  R 
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3.  The  Steklov  problem 


Coming  back  to  (2.9)  we  introduce  the  Steklov  eigenvalue  problem:  Find 
all  pairs  (S.A),  SeH^CSlj^)  such  that  (2.9)  holds. 

In  the  case  that  =  0,  the  trivial  function  S  =  1,  (A  =  0)  will  not 

be  considered  here. 

We  C2U1  cast  the  Steklov  problem  in  a  different  form.  To  this  end 
define  the  operator  T:  — »Hp(£2j^)  such  that  (with  (2.10)) 

(3.1)  Bjj^(Tu,v)  =  bj^(u,v)  Vu.veH^dlpj) 


The  operator  T  is  obviously  selfadjoint  and  is  well  defined  by  the 


1  2  0 

coercivity  of  .  Since  the  trace  mapping  u  i - >  u|j.o,  H  (0^^)  - >L  (T^^). 

R 

2  0  2  0 

is  compact  and  the  bilinear  form  bj^(u,v)  is  continuous  on  L  (F^)  xL  (F  ), 

the  operator  T  :  compact. 

Now  the  Steklov  problem  (2.9)  can  be  cast  in  the  following  form: 

1, 


Find  (S,A),  9  6  11^(0^^)  such  that 

(3.2)  ATS  =  S 


or  denoting  A  =  t 
A 


(3.2a) 


TS  =  AS. 


It  is  known  from  the  theory  of  compact  selfadjoint  linear  operators  that 


there  au-e  countably  many  eigenpairs  (S.,A.),  A  .  e  R,  with  no  accumulation 

0  3  0 

points  except  at  A  =  0.  Furthermore  the  eigenfunctions  S.  yield  an 

0 

orthonormal  basis  of  the  closure  of  the  range  of  T.  If  A  is  a  simple 

\J 

eigenvalue,  then  S^  is  uniquely  defined  up  to  a  multiplicative  factor.  For 

multiple  eigenvalues  A  =...A.  we  will  assume  that  S,,---,S,  au'e 

J  J  J+m 

orthonormal  with  respect  to  B_  .  In  this  case  only  the  span  {S,,---,S.  }  is 

Wr  j  j+m 
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unique  but 

the  orthogonal  basis  (S.,*-' 

\J 

is  not  unique 

From 

(3.1)  we  readily  see  that  for 

u 

=  Tv  we  have 

(3.3a) 

Au  =  0 

in 

(3.3b) 

u  =  0 

on 

^R,D 

(3.3c) 

on 

^R,  N 

(3.3d) 

du 
an  ~ 

on 

Hence  the  eigenfunct: 


(defined  in  the  variational  sense); 


£(£2^^)  =  {ueHpCnp^),  Au  =  0  in  £2^^,  u  =  0  on 


on 


^R.N>- 


We  also  see  that  functions 


S.|_o  form  an  orthogonal  basis  of  the  space 
R 


As  usually  we  will  assume  that 


(3.4a) 

sold 


B„  (S..S,  ) 

Hr  J  k 


'  0 

1 

V 


for  k  *  j 
for  k  =  j 


(3.4b) 


0  for  k  *  j 

i  =  A  for  k  =  J 

Aj  J 


In  the  case  that 

K 

modulo  constants. 

So  far  we  assumed  that  Hp(£2j^)  is  a  real  space.  We  will  extend  it  to 
the  complex  space  C  by  defining,  as  usual. 


2  0 

=  <f>,  we  will  understand  ^(£2^,)  and  L  (r„)  as  spaces 
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(\7u)  •  (7v)dx, 


-1 


with  b^(u,v)  being  analogously  defined. 

Denote  by  p(T)  the  resolvent  set  of  T,  i.e.  p(T)  =  iz | z € C, (zI-T) 
exists  as  a  bounded  operator  on  (by  I  we  denoted  the  identity 

operator).  Further  let  <r(T)  be  the  spectrum  of  T  i.e.  trCT)  =  C-p(T). 

For  any  zep(T)  denote  “  (zI-T)  the  resolvent  operator. 

Let  p  be  a  nonzero  eigenvalue  of  T  with  multiplicity  m+1  and  y 
be  a  circle  in  C  centered  at  p  which  lies  in  p(T)  zmd  which  encloses  no 
other  point  of  (r(T)  than  p.  Then  the  spectral  projection  associated  with 
T  at  p  is  defined 


(3.5) 


E  =  ^ 


R  (T)dz 
z 


see  e.g.  [BOl].  Now  we  prove 


Theorem  3.  1  Let  (Sj^.A^^)  be  a  Steklov  eigenpair.  Then  Sj^  has  the  form 


(3.6) 


Sj^(r,0,v)  =  r  Wj^(0,^) 


Proof  Let  Sj^(r,0,v)  be  the  Steklov  eigenfunction.  Then  we  can,  for  any 
0  <  r  ^  R,  define 


(3.7) 


S  (r,0,  ^)S  ,(R,  0,  (ip)ds 

K  J 


R 

For  0  <  r,  <R  apply  now  Green’s  formula  for  S,  (r,0,v)  aind  SJ — r,e,<p) 
1  K  j  r^ 


in 


Then  it  can  be  readily  seen  that 
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(3.8) 


,rR 


B  (Sj^Cr.a.v).  s 

«!  1 


r  as. 


^  ir.9.<p).  Sj(R,e.v>)ds. 


fO 

ri 


On  the  other  hand  we  have  by  scaling  (2.9) 


(3.9) 


,rR 


1 


Further  we  see  that 

(3. 10) 


=  A.  r~^ 
J  1 


Sk  (r^.e.^p),  Sj(R,e,(p)ds 


r° 

ri 


(k),  w  -1 

=  a.  (PjjAjr,  . 


r  as. 


^  (r^.e.^)).  s^(R,e,»))ds. 


fO 

ri 


=  ^ _ (r  ) 

dr  ^  r 

where  the  derivative  on  the  right  hand  side  of  (3. 10)  is  understood  in  the 
weaJc  (distributional)  sense.  Combining  now  (3.7)-(3.10)  we  get 

.  (k) 

J  -  A  r-' 

dr  J 


and  hence 
(3. 11) 


a^^^(r)  =  0<r<R 

\)  0 


(k) 

From  (3.4b)  we  get  =0  for  k«J. 


Since  S,|_o  is  a  basis  of  L  (F^)  we  obtain 
J  Tr  R 
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1  \  -  \ 
Sj^(r,0,^)  =  —  -  Sj^(R,0,v)  =:  r  Wj^(0,y) 

R 


which  was  to  be  proven. 
We  get  Immediately 


Corollary  3. 2.  The  singularity  function  au'e  exactly  the  Steklov  eigen¬ 
function  Sj  (up  to  a  factor  or  lineair  combinations  for  multiple  eigenvalues). 
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4. 


The  vertex  intensity  factors 

Let  u  be  the  solution  of  the  boundary  value  problem  (2.1)  (or 


equivalently  (2.5))  ajid  assume  that  (2.6)  holds.  Then  we  have  u| €  i?(£lj^) 
and  hence  by  corollary  3.2  u  can  be  written  in  the  terms  of  the  basis 
functions  (Sj,S2,...> 


(4.  1) 


J=1 


Let  us  show  that  (4.1)  implies  the  decompxjs i t i on  (2.11)  with  C  =  C.. 

U  vJ 

Let 
(4.2) 

and  for  0  <  p  ^  R 


u_  =  u  -  y  C.S.  =  y  C.S. 
0  Z-  J  J  J  J 


A  +-<s 

J  2 


A  .+->s 
J  2 


F(p):  = 


|7Uq)  dx. 


a 


then  by  the  orthonormality  of  S .  we  get 

\) 


F(R) 


=  ^c5  =  c<.. 


^*->S 

J  2 


Further  let 


“o'*’  '  “o(r*)  =  Z  ^j(r) 


A.+->s 
J  2 


then  for  sufficiently  small  e>0 
(4.3) 


2Aj  f  ,1+2(s-1/2)+e 


<  e 

■  R 


z 


x2  .  -  2s+e 
Cj  <  c  p 


Aj*5>s 


A  +->s 

J  2 
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Let 


■II 


-2s  2  2 

r  “  I  7Uq  I  ^rdr  du> 


0a,erO 

where  d«  denotes  the  surface  element  on  F?.  Hence  integrating  by  pairts 

K 


we  get 


-2  -2s  -2-2s 

G  =  R  r  '^^F'Cridr  =  R  "^®F(R) 


-R 

-2  -2s- 1 

+  2sR  r  F(r) 


and  using  (4.3)  yields 


-2s  2  -2  -2-2s  -2  -2s- 1 

r  |7Uq|  dx  =  R  G  =  R  F(R)  2sR  r  F(r)dr 


5  C(1  +  2s  r'^'^^^dr)  S  C. 


Therefore  (2.11b)  holds  with  Cj  =  C^.  Hence  we  can  use  (4.1)  and  (3.4)  to 
express  the  intensity  factors  C, 


(4.4) 


-j  =  ^  f  “ 

J  nO 


Replacing  u,Sj,Aj  by  their  finite  element  approximations,  we  will  use  (4.4) 
in  section  7  to  compute  Cj  numerically. 


Remark  4.  1  Note  that  in  the  ceise  when  is  a  multiple  eigenvalue  the 

functions  S.  are  not  unique  and  hence  C.  depends  on  the  choice  of  the 

singularity  function  S  .  For  a  simple  eigenvalue  A.  the  eigenfunction  S. 

J  3  3 

is  unique  up  to  a  factor  of  -1. 
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5. 


Finite  Element  method 


Let  ?*q(n),q  =  1,2,  •••  be  the  partition  of  (2  into  the  set  of  open 

elements  i  =  1,2,  •••  M(q),  be  the  aissociated  finite  element 

C--spaces  and  let  N(q)  be  the  dimension  of  W  (12).  Further  let  X  = 

O  q  q 

{9  (12), W  (0))  and  M  =  {X  }. 

q  q  q 

Denoting  by  the  finite  element  solution  of  the  problem  (2.1) 

satisfying 

(5.1)  B„(u  ,v)  =  B_(u,v)  Vv€W  (£2) 

£2  q  £2  q 

we  obviously  have 


(5.2) 


u-u 


q 


where  B(u,v)  and  ||  •  were  defined  in  Section  2  (see  (2.3),  (2.4)).  If 

the  data  8.g|^  in  (2.1)  are  sufficiently  smooth  say  geH  ^*®(£2),  gj^  = 

for  some  GeH^  “(£2),  s>0  then  the  weaik  solution  H€Hp(£2)  of  (2.1) 

1  +s  ^  ^ 

(resp.  (2.5))  exists  and  is  unique.  It  belongs  to  H  (£2)  where  £2  = 

£2-  U  E^,  being  the  A  neighborhood  of  E..  In  the  neighborhood  of 
j=l  J  J  J 

Ej,J  =  l,-'',n  and  =  0,l,**-,m,  the  behavior  of  u  is  determined  by 

the  edge  auid  the  vertex  singularity  functions.  The  approximabl 1 ity  of  these 
functions  determines  the  convergence  rate  of  finite  element 

solution.  We  assume  that  the  convergence  of  the  finite  element  solution  is 
characterized  by  a  nonincreasing  function  F  :  IN — >IR^,  with  F(q)  — >0  as 

q >00.  More  precisely,  we  will  say  that  u^  is  F-convergent  if  there  exists 

a  constant  C  independent  of  q  but  depending  on  u  such  that 


(5.3) 


-  CF(N(q)) 

^  <en^(n)  ^  ^ 
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Remark  S.l:  Later  we  will  aissume  that  an  estimate  of  the  form  (5.3)  with  the 


same  function  F  also  holds  for  a  class  of  solutions  which  will  be  specified 
in  Assumptions  A^,  A^,  A^  below.  It  may  be  that  the  function  F  which 
satisfies  Assumptions  A^,  Ag,  A^  gives  a  less  them  optimal  error  estimate  in 

(5.3) ,  e.g. ,  in  the  exceptional  case  when  the  solution  u  is  smooth. 

An  interesting  question  is  the  characterization  of  all  functions  u 
satisfying  (5.3)  for  a  given  function  F(N(q))  eind  a  given  sequence  of  meshes. 
For  the  description  of  such  class  of  functions  in  a  particular  case  we  refer 
to  [BKP].  In  practice  we  can  assume  that  “  CF(N(q))  which  is  a 

typical  case,  but  in  the  sequel  we  will  only  assume  that  (5.3)  holds. 

Let  us  describe  the  convergence  function  F  for  some  typical  examples. 

Example  5.  1 .  The  h-version  method  on  a  quaisiuniform  mesh.  Let  t)e  the 

standard  family  of  quaisiuniform  simplical  (in  general,  curved)  meshes  of  the 
size  ^  (see  e.g.  (C1],(C2]).  Let  W^dij^)  be  the  set  of  functions 
which  are  polynomials,  of  total  degree  <  d  on  each  simplex.  Then  we  have 

(5.4)  F(N(q))  =  N(q)"^^^ 
where 

(5.5)  3  =  min(d,  s,  <r-G} ,  g>0  arbitrary 

(5. 6)  <r  =  min-^A|*^^  +  =  0,---,m,  t  =  1,2, •••,n 

(k) 

where  A^  is  the  smallest  vertex  singularly  exponent  for  the  vertex 

singularity  function  (see  (2.11))  and  where  is  the 

u 

internal  argle  of  Q  at  the  edge  F^.  If  E^  =  F^nFj  auid  i  €  D,  JeN 

or  i  €  N,  JeD  (i.e.  the  Dirichlet  condition  is  prescribed  on  one  side  and 

Neumam  condition  on  the  other  side  of  E^)  then  instead.  By 

2u 
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1  ^0*”G 

the  regularity  theory,  we  have  then  ueH  (£1),  e<0  arbitrsiry  where  we 
1 

denoted  by  H  (0)  the  st2utdard  Sobolev  space  with  fractional  derivatives 
(for  definition  see  e.g  [BL]).  (5.4)  then  follows  from  the  standsurd  theory  of 
the  finite  element  method. 


Example  5.2.  The  p-version  of  the  finite  element  method.  Here  is 

fixed  mesh  of  simplices  (generally  curved)  zund  W^fO)  Is  the  space  of  all 
functions  in  H^(il)  which  are  polynomials  of  degree  q  on  each  simplex. 


Then  we  have 

(5.7) 


F(N(q))  =  N  g)/3^  e>0  arbitrary 


where 

(5.8) 


8  =  min(s,2<r), 


where  <r  is  given  in  (5.6)  (see  [Dl],  ID2]). 


Exzunple  5.3.  The  h-p  version  of  the  finite  element  method.  Here  P  (R)  is 

-  q 

sequence  of  properly  selected  meshes  zuid  is  the  space  of  functions  in 

H^(R)  which  aire  polynomials  of  degree  R(q)  with  R(q) — >oo  as  q — 
properly  selected.  We  assume  that  g  is  an  ainalytic  function  on  R  and  g 
is  an  analytic  function  on  every  face  r^.  We  ceui  then  expect  that 


N 


(5.9) 


F(N(q))  =  e"^^ 


1/4 


(5.9)  was  proven  in  the  caise  that  RcR^,  n  =  1,2  when 

F(N(q))  =  e 

see  [BC],  (GP). 

We  need  now  to  make  additional  assumptions  about  the  family  M. 

Ai:  Fp  coincides  with  the  boundaries  of  the  elements  of  the  partition 

P  (R). 

q 
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Let  u€Hp(£l),  Au  =  0  on  and  Au  =  0  on 
u.  =  u|_  2md  u_  =  u|  assume  that  u,  resp.  u„  can  b« 

extended  to  a  harmonic  function  on  resp.  Q.-ZL^^  and 

on  Fd^  Then 


Denoting 

(auialytically) 

=  0,  t  =  1,2 


(5.  10) 


inf  llu-CII  i  CF(N(q)) 
<€W^(n) 


Let  ueifCn^i^)  (as  defined  in  Section  3).  Then 


(5. 11) 


inf  llu-<ll_.„  .  <  CF(N(q)). 
<€W^(nj^) 


In  (5.10)  auid  (5.11)  the  function  F(N(q))  is  aissumed  to  be  the  same  as  in 
(5.3). 

Let  us  now  discuss  the  validity  of  the  assumptions  A^,  A^,  A^- 

i)  Assumption  Aj:  It  can  be  satisfied  by  the  stamdard  finite  element 
technique  using  the  binding  mapping  of  the  maister  element  (for  the  h-p  version 

see  [BG]  for  details). 


ii)  Assumption  A3:  Here  u€if(n2f^)  implies  that  the  function  u  has  only  a 
vertex  singularity  at  =  0  aind  the  edge  singularities  at  the  edges 
Ej^,  i=l,  ■  • ' , n^,  but  it  hais  no  additional  singularities  at  Therefore,  the 

assumption  A3  is  satisfied  in  the  examples  5.1,  5.2  and  5.3  by  the  direct 
application  of  the  corresponding  approximation  results. 


ili)  Assumption  A2:  We  will  briefly  sketch  its  validity  for  the  examples 
mentioned  above. 


Example  5. 1.  The  h-version  on  the  uniform  mesh.  For  simplicity  we  will 
restrict  ourself  to  the  case  d  =  1.  Consider  the  space 
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«2  = 

Ilull^  =  ||u||\  +  llull^ 

Ha  ircnR)  ir(n-nR) 

Further  let  be  the  energy  projection  operator  of  on  W^(n).  Then  by 

the  standard  approach,  see  e.g  [Cl],  (C2],  we  have  for  v e 

(5.12)  '''-V'E(a)^™l'''H,(n)- 

Applying  the  interpolation  theory  (see  e.g.  [BL])  we  see  that  for  1 < s < 2, 

0  =  s-1 

and  hence  for  any  u  e  , 

llu-nhUii^(„).ch®lluii^. 

Assume  now  that  u  satisfies  assumption  Aa-  Then  u  has  the  same  type  of 
singularities  in  12^^  and 
Hence  u| with  «r  as  in  (5.6)  and 
(5. 10)  follows. 

Example  5.2.  The  p-version.  Here  the  validity  of  A^  follows  by  applying 
Dorr’s  results.  (see  [Dl],  [D2]).  He  approximates  first  the  solution  u 
having  singulairity  of  the  vertex  and  edge  type  element  by  element  using 
weighted  (Legendre  type)  spaces,  imposing  continuity  at  the  vertices  of 
elements.  Then  the  difference  of  the  approximation  on  the  edges  and  faces  of 
the  neighboring  elements  (the  discrepancy)  in  the  (Legendre)  weighted  spaces 
norm  is  estimated.  Then  it  is  shown  that  it  is  possible  to  extend  this 
discrepancy  into  the  elements  where  the  extension  mapping  is  continuous  from 


as  solutions  of  (2.1)  with  smooth  g  and  g^^. 
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the  weighted  spaces  on  the  boundary  to  in  the  element.  Hence  only  the 
smoothness  of  the  solution  in  each  element  is  employed.  Realizing  that  in 
every  element,  the  function  u  can  be  decomposed  into  a  smooth  function  and 
the  singuleu'  functions  the  arguments  of  [Dl]  and  [D2]  apply. 

2 

Example  5.3  The  h-p  version  in  IR  .  Here  the  singularity  is  only  in  the 
vertices  of  the  domain  and  there  are  only  finite  number  of  elements  which  have 
boundary  on  F^.  As  in  [BG]  we  approximate  u  separately  on  every  element 
and  then  remove  the  discrepancy  (discontinuity)  of  the  approximation  on  the 
boundary  of  elements.  Because  the  solution  in  every  (closed)  element  which 
nonempty  intersection  with  F^  is  analytic,  the  arguments  used  in  [BG]  are 
immediately  applicable. 
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6.  Computation  of  the  Singularity  Functions 

The  singularity  functions  Sj  (see  Corollary  3.2)  will  be  computed  by  the 
finite  element  method  as  the  approximate  eigenfunctions  of  the  Steklov 
eigenvalue  problem.  We  will  assume  that  the  assumptions  A  -A  introduced  in 

X  O 

the  section  5  hold. 

The  finite  element  solution  of  the  Steklov  problem  is  based  on  the 
variational  formulation  (2.9)  (2.10):  Find  €  W^(nj^)  and 

such  that 

(6.1)  Bjj^(S^‘’^v)  =  Aj‘’^bj^(S^‘^^v),  V  veW^(np). 


By  normalizing  the  eigenfunctions  and  orthogonal izing  the  eigenfunctions  for 
multiple  eigenvalues,  we  have,  analogous  to  (3.4  ab). 


(6.2a) 


(6.2b) 


-tql 

•^k 


)  = 


for 

for 


for 

for 


k*  J 

k  =  j 

k*  J 
k  =  J 


Let  further  (Sj,Aj)  J  =  1,2, ••• 
(sec  [BOl]). 


be  the  exact  eigenpairs. 


Then  we  have 


Theorem  6.1.  Let  S.,S.  ,,**-,S.  be  the  eigenfunctions  associated 

J  j+l  J+n> 

eigenvalue  A,  with  wu.lt  iplicity  m+1  (i.e.  A=A.  ,  =  ...A.  ). 

J  K  J  J+l  J+m 

q  sufficient iy  large  there  exist  exact  eigenfunction  S.  .  ,l=0,-- 

J+»-»  q 

(depending  on  q)  and  satisfying  (3.4)  such  that 


(6.3) 


'^J+«.q  ^j+f"E(nj^)  - 


(6.4) 


|A.  <  Cc^ 

J+«  J+<  J 


to  the 
Then  for 
•  m 
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where 


(6.5) 


e ,  =  sup  inf  IIS  -  Cllpfo  ^ 

<€W^(«r) 


and  M(A.)  is  the  eigenspace  associated  to  the  eigenvalue  A.  of  multiplicity 
J  J 

m  +  1.  o 

We  will  show  later  that  we  can  impose  additional  conditions  on  S. 

J  +  C,  cj 

The  constzmt  C  in  (6.3)  (6.4)  is  independent  of  q  but  depends  on  various 

other  factors  (see  (B02]  for  the  discussion).  In  the  sequel  we  will  also 

write  S.  instead  of  S .  if  no  misunderstanding  occurs. 

J  J.  q 

With  the  assumption  we  get 


(6.6) 


e  .  <  C.F(N(q)) 
J 


and  hence  from  the  theorem  6.1  we  get 


(6.7) 


(6.8) 


lA.-A^'^^l  <  C(F(N(q)))^ 


Remark  6.1  Note  that  the  bilinear  form  bj^(u,v)  in  (2.10)  depends  only  on  u 

and  V  on  F?-  This  allows  us  in  practice  to  eliminate  first  all  unknowns 

R 

inside  aind  obtain  an  eigenvalue  problem  on  F^  only.  By  proper  orderings 
of  unknowns  the  LU  decomposition  of  the  stiffness  matrix  on  £2^^  can  be  used 

in  the  computation  of  u^  i.e.  for  the  computations  of  the  finite  element 

solution  of  (2.1). 

Let  us  now  write  the  finite  element  solution  of  the  Steklov  eigenvalue 
problem  in  a  different  way,  which  is  the  basis  of  the  estimates  (6.3)  (6.4). 


This  will  be  used  in  the  next  section  too. 


Let  be  the  elliptic  projection  of  H^dlj^)  into  W  (£1^^)  defined 
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(6.9) 


Define  now  T  =  where  T  is  given  in  (3.1).  Then  T  converges  to  T 

q  Q  q 

in  the  norm  of  linear  mappings  11^(0^^) — >H^(nj^)  as  q — »«. 

Further  consider  the  spectral  projection  E(Aj^)  onto  the  eigenspace  of 


A,  of  T  given  by  (3.5)  where  y  is  a  circle  which  encloses  A, 
J  J 


but  no 


A^^  which  Aj^^A.; 


(6. 10) 


E(A“h  =  ^  (T-zD  ^dz 


Let  S.  ^  be  an  eigenfunction  of  (6.1).  This  implies  by  definition  of  T 
J  ‘ 

that  is  an  eigenfunction  of  T  : 

J  ^ 

q  J  j  J 

We  have  then  the  following  relation  for  the  projection  of  an  eigenfunction 


of  T^  onto  the  eigenspace  of  T 


(6. 11) 


-E(Abs^.^^  =  ^  (T-zI)  ^(T  -T)(T  -zl ) "^S^.^^^dz 
j  J  J  2ni  q  q  J 


In  (6.8)  the  roles  of  T  and  T^  are  reversed  in  compEU'ison  with  [BOl], 


section  7. 


Let  us  impose  additional  conditions  on  S .  : 

J.  q 


Theorem  6.2.  Let  I  =  0,...m(j)  be  the  approximate  (finite  element) 

eigenfunctions  associated  with  the  eigenvalue  A^  of  multiplicity  m(j)  +  1. 

Then  there  exists  S,  .  eM(A  )  satisfying  (3.4)  such  that 

j+c,  q  J 


(6. 12) 


'^J+«.q'^J+^’E(DR)  -^*'j 
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(6. 13) 


IS.,  -E(a"^)S^,^1||_,„  , 
j+i. q  J  j+i  ecar)  j 


where  is  defined  by  (6.5). 


Proof.  As  mentioned  previously,  E(Aj  orthogonal  projection  of 

[q] 


into  M(Aj).  Let 


(6. 14) 

We  have  by  orthogonality 


(6. 15) 
and  hence 

(6. 16) 


lglq]||2  _  ..g  ..2  ..g  _  „[q]..2 

'^J+«"E(nR)  "^j+£.q"E(QR)  "  J+i.q  ^j+i"E(nR) 


Now  we  will  construct  by  induction  the  orthonormal  system  of  functions 

S,  -  =  0, satisfying  (6.12)  (6.13)  by  the  Gram-Schmidt 

j+c,  q 

orthogonal izat ion  process.  Assume  now  that  (6.12)  and  (6.13)  hold  for 
^  =  0,  •••,  k-1  (for  k  =  0  no  assumption  is  necessary  since  S, 

S 

j+K,q 

k-1 

^j+k,q  ■  ^J+k,q  q’ ^J+£,  q]^j+«,  ^ 


>k.q 

(6. 17) 


J+k,q 

auid  only  eq.  (6.18),  (6.21)  below  Eire  used).  Then  we  define 


1=0 


and 
(6. 18) 

Consider  first 


S 

J+k.q 


J+k.q 


>k,q‘'E(nR) 


k-1 


^J+k,q  ~  ^J+k,q  X  ®S2R[^j+k.q’^J+^.q]^j+^.' 

1=0 


Noting  that  by  (6.14),  (6.2a)  we  have 
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B, 


I  fs  -  s 


><.q]  " 


k,^  =  0. 


and  hence 


Therefore 


This  yields 


®nR 

ctq]  cfqil  - 
[^j+k-^j+£j  ■ 

0.  k.^ 

=  0.  • 

"  *  » 

m(  J) . 

-  s 

§  ,  1 

=  0. 

k.  1 

=  0, 

^R 

[  J+k  J+k.q 

J+«.  qj 

^  » 

% 

• 

_^J+k.q’^J+£.q 

]  [^J+k.q’^ 

.[q] 

] 

=  B 

fs  ^ 

np 

L  J+k.q  ""j+k 

’  J+«J 

=  B 

s  . 

1  . 

Hr 

L  J+k.q  '’j+k 

*  j+^ 

j+^.q 

J 

I  Bo  |s,^.  ,  hi  Cc^. 

ArL  J+k,q  j+£,qJ  J 


'j+k.q"  Vk,q"E(IlR)  * 


and  by  induction  hypotheses  we  get  also 

(6.20) 


^J+k,q  ^j+k,q'*E(nj^) 


Using  (6. 16) ,  we  get 

(6.21) 

and  from  (6.20)  we  get 
(6.22) 


Hence  from  (6.18)  and  (6.20)  we  get 


'^J+k,q  Vk.q"E(nR) 


which  was  to  be  proved. 


,  m(  J) 

k  *  < 

• • , m( j) 
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7.  Computation  of  the  Vertex  Intensity  Factors 

The  exact  value  of  the  vertex  intensity  factor  is  given  by  (4.4). 
we  will  numerically  compute  as  follows: 


(7. 1) 


We  then  want  to  estimate  the  error 


u 

q  J 


ds. 


Hence 


As  we  said  in  the  remark  4.1,  the  stress  intensity  factors  Cj  depend  on 
the  choice  of  exact  eigenfunctions  Sj.  The  eigenfunctions  are  not  uniquely 
defined  if  the  associated  eigenvalue  A.  is  not  simple.  In  this  case  we  will 
consider  the  stress  intensity  factor  C.  .  which  are  associated  with  the 
eigenfunctions  ^j+£  q’  ^  ~  0,*"'m,  which  were  introduced  in  theorem  6.2.  The 
decomposition  (2.11a)  of  the  exact  solution  can  be  expressed  in  this  baisis 
with  vertex  intensity  factors  C,  depending  on  q: 

j»  q 


u  =  u„  +  7  C ,  § , 

0  4.  J.q  J.' 


A>-<s 

J  2 

For  simplicity  we  will  write  in  what  follows  S,  instead  S. 

J  J.q 

We  have  now 


Theorem  7.1  Assume  that  the  assumptions  A^^  -  A^  are  satisfied.  Then  for 
sufficiently  large  q  we  have 

(7.2)  -C  I  <  C(F(N(q)))^ 

J  J  »  ^ 

Proof.  Using  (7.1)  we  have  to  show 

(7.3)  |B„  fu  -B„  (u.S,)|  <  C(F(N(q)))^ 

q  J  J  Wr  j 

We  h, 
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We  will  estimate  D^.i  =  1,2,3  separately. 

1)  Estimate  of  D^.  Using  assumption  A^,  we  get 


(7.5) 


ID3l^||Sj-s''’'l|,,^,l|u-u^,^,„^,  S  CF(N(q))2 


which  follows  immediately  from  (6.12). 


2)  Estimate  of  D^.  Let  weHjjCQ)  satisfy 


(7.6) 


B-(w,v)  =  (S.,v)  Vv€H;.(n) 

n  wr  j  D 


The  function  w  exists  auid  is  unique.  It  can  be  readily  seen  that  (7.6) 
is  the  variational  formulation  of  the  problem 


(7.7) 


Aw  =  0 

in 

«R 

Aw  =  0 

in 

wt  =0 

D 

dw 

1  =  ° 
‘'n 

(w]  o  =  0 

R 

rdwi 

L^J 

_  a 

_o  dn 

R 


By  [*1-0  we  denoted  the  jump  across  f^.  Let  w 


“2  '  “'n-ii, 


define  w  in  n„„  by 
1 


(7.8) 


“l  = 


w. 


W  ”  — — 

2  1+2A 


in 

:h  -  “2b-Sr' 


Because  we  have  S,  =  on  P?,  we  have  [w  ]  o  =  0  and  since 

0  0  ^  ^ 
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a  L  d1+2Aj„  1  2A.+1  as.  *  fs  ~  1 

^  S ,  -  R  ■’S  ,  =  — i —  1-  we  get  w.  „  =  0 

an[^  J  -Jj  Aj  an  ®  [an  ij  0 

R 


Further,  Aw^  =  0  on  and  AW2  =  0  on  02^^  - 


Hence 


AWj  =0  on  shows  that  w  can  be  analytically 

extended  from  into  n2f{  a  harmonic  function. 

Analogously,  we  can  extend  the  function  W2  on  by  defining 

'  ”2  in  n  -  Qj^ 


(7.9) 


“2  = 


“1  *  “b  -  ^2 


Because  of  (5.1)  we  have  for  any  w  eW  (Q) 

<1  q 


D„  =  (S,,u-u  )  =  B(w,  u-u  )  =  B(w-w  ,  u-u  ), 

2  Qr  J  q  ’  q  q  q 


Using  assumption  A2  we  get 


(7.10)  ">2l  ^CF(N(q))‘ 

3)  Estimate  of  D, 


<€Wj(n) 


Define  as  in  (6. 14) 

(7.  11)  =  E(a'^)sJ‘^^ 


Then 


J  J  ■  J 


For  sufficiently  lairge  q  we  have  using  (6.11) 


-  E(a"MSj^^ 


1 

2irl 


(T-zI)"^(T  -T)(T  -zD'^S^'^^dz. 

q  q  J 


Since  T 

q  J 


's''’'  we  have 
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and  hence 


q  J  ,  J 


(A-M  ^-z 


'  2iti  I  ®n»(' 


(T-2l)  '(T  -T)  — r-T — i —  s''’'.uldz 


By  (3.1)  we  have 


and  hence 


which  yields 


Vw.veH^CQj^) 


B-  ((T-zI)'^w.v)  =  B_  (w. (T-zI)  ^v) 
“R  Hr 


Because  u|_  eifCno)  we  can  write 
*1r 


(7. 12) 


(7. 13) 


Hence 


Let  us  define 


(7. 14) 


ic,r<». 


-zi)-'u  =  £(a;;^-z)-'c^s^ 


^  ^  -^]  (T-z)"^u  dz 

'  I  ((''J'’')  H  (\’-^r'dz. 

k=l  ^7 
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Since  ~  ~  is  inside  of  the  circle  jr  and  all  other 


are  outside  we  have  from  the  residual  theorem  using  the  poles  at 
z  =  auid  z  =  ^ 


Therefore 


0  for  k  =  J+^,  t  =  0,  ■  ■  • ,  m(  j ) 

^  k  j+£,  <  =  0,---,m(J) 


(7. 15) 


“  =  Z  ‘c^Sj 


k=l 

k*j+« 


and  hence  for  any 


(7. 16) 


e'^ECnR)  *  "^~^''E(nR)' 


Let  us  estimate  both  terms.  To  this  end  we  write 

=  ‘V’V* 


(T^-T)S.  ,  =  (n^-I)TS.^,  =  (II^-I)a'^S.^,. 
q  J-*-^  q  j+«  q  j  j+« 


Then  using  assumption  A^  we  get 


c=u,  •  •  • ,  ml  J  J 


Further  because  W^q~^llE({j  )-)E(n  )  ”  ^  (for  q  sufficiently  lairge)  we  get 
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(7. 17) 
aaid  hence 

(7. 18) 


Let  US  now  estimate  the  second  term  In  (7.16). 

We  know  that  uei?(f22j^).  Hence  u  (x)  =  u(2x)  €  JPCHr) .  Therefore  we 


can  write 

(7. 19) 
with 

(7.20) 

On  the  other  hand  we  have  on  OR/a- 


-z 


c^s. 

k  k 


k=l 


1^  <  00. 


k=l 


(7.21) 

and  (7.19)  yields 

(7.22) 

(7.23) 


u  (x)  *  u(2x) 


=  ^C,S,(2x)  =  y  C.Z^'^SJx] 


k  k' 


k=l 


C*  =  C.  2^“ 
k  k 


/L.  k  k' 

k=l 


k=l 


(7.24) 


Now  we  consider  u  (x)  =  u(2x).  We  get  from  (7.15)  for  x  e  Qr 

S-(x)  =  f , 


c»,,S.  (x) 
k  k 


k=l 


where 


(7.25) 


J 

0 


k  =  j+^  . 
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Since  for  k  *  J+i  the  numbers  8U'e  outside  of  y  and  ^  is 

inside  of  y  (for  q  sufficiently  large)  we  have 


(7.26) 


(7.25).  (7.26),  (7.22),  (7.23)  yield 

k=l  0  k=l 

Therefore  u  and  hence  u€if(n2f^)-  Using  assumption  we 


have 

(7.27) 


inf  liu-CII  <  CF(N(q)). 

<€W^(nR) 


Hence  from  (7.16)  usin<i  (7.18)  and  (7.27) 


(7.28) 


Using  now  (6  13)  we  get 


<  CF(N(q))‘ 


i.e.  IDjI  S  CF(N(q))  what  was  to  be  proven. 


Remark  7.1  We  have  assumed  that  (2.6)  holds.  This  allows  us  to  use  the 
decomposition  (2.11).  In  the  general  case  the  decomposition  is  more  complex. 
The  present  theory  can  be  generalized  to  this  ceise  but  we  will  not  address  the 
problem  here. 
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8.  Numerical  Results 

Let  us  consider  the  following  mixed  boundary  value  problem.  Let  £1  be  a 
cube  with  a  crack, 

Q  =  {|Xj|  <  1.  J  =  1.2.3>  -  {(x^.X2.0)|0<x^  <  1,  O^x^^  1> 
as  shown  in  Fig. 8.1. 


Figure  8. 1.  The  cracked  cube 

On  both  sides  of  the  crack  auid  the  two  faces  x^  =  1  and  x^  =  1 
Dirichlet  conditions  u  =  0  aire  prescribed.  On  the  remaining  4  faces  of  the 
cube  the  following  Neumaain  conditions  aire  prescribed 

i)  =  0  on  the  faces  x,  =  -1,  x„  =  -1 

on  1  2 

ii)  ^  =  cos  ^(x,  +  l)  cos  ^(x.,+  1)  for  x  =  ±  1 

on  4  1  4  2  3 

The  boundsu'y  conditions  imply  that  the  solution  will  not  have  singulari¬ 
ties  at  the  edges  and  vertices  of  the  cube  (since  one  can  use  even  and  odd 
extensions).  The  solution  has  a  singularity  at  the  two  edges  of  the  cracks, 
(0, 1 )  X (0)  X {0>  and  {0}  x  (0, 1 )  x {0>  ard  at  the  vertex  at  (0, 0,0).  For 
a  more  detailed  description  of  the  edge  amd  vertex  singularities  in  this 
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example  see  [PSl],  The  (smallest)  leading  eoge  singularity  exp>onent  is  a  =  |. 
The  leading  vertex  singularity  exponent  is  not  analytically  knovm,  our 

computational  results  below  show  »  0.2966.  We  use  symmetry  and  analyze 
only  the  half  cube  (-1<X2<0)  computationally.  The  p-version  of  the 
finite  element  method  implemented  in  the  program  STRIPE  [S]  is  used.  We  used 
two  different  meshes,  the  "unrefined  mesh"  and  the  "refined  mesh. "  Both 
meshes  contain  a  ball  around  the  vertex  at  0. 

Figure  8.2  shows  the  unrefined  mesh.  For  clarity,  the  mesh  inside  is 

shown  separately. 


Figvire  8.2.  The  basic  unrefined  mesh  around  the  vertex  at  0. 

In  Figure  8.3  we  show  the  refined  (geometrical)  mesh  with  6  layers. 

The  sizes  of  adjoining  layers  have  a  ratio  of  <r  =  0. 15  for  the  mesh  used  in 
the  computation  (in  Fig.  8. 3  o"  =  0.5  is  used  for  clarity). 


Figxire  8.3. 


Detail  of  the  geometrically  refined  mesh  in  the  neighborhood 
of  the  vertex  at  0. 
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and 


In  table  8.1  we  repxjrt  N(q)  and  the  computed  values  ol' 

2 

||u  llr-rnA  the  mesh  shown  in  Figure  8.2  as  function  of  the  degree  p  of  the 

q  t. I  li  ) 

elements.  Table  8.2  gives  analogous  results  for  the  refined  mesh  (with  6 

layers)  shown  in  Figure  8.2.  From  tables  8.1  sind  8.2  we  cleau'ly  see  the  effect 

of  the  refinement  of  the  mesh.  By  using  a  more  refined  mesh,  a  higher  value 

2 

of  p  and  extrapolation,  we  estimate  the  exact  values  for 


(8.1)  =  0.29658  =  10.123,  =  8.66908 


Table 

2 

8.1.  The  values  of  ^^>^2  ^ind 

the  mesh  shown  in  Fig 

P 

N 

^1 

ptq] 

^1 

''^"e(£2) 

2 

133 

0. 315743 

9,50571 

8. 5697572 

3 

231 

0.313768 

9.55731 

8.5867675 

4 

420 

0.307749 

9.75874 

8.6263436 

5 

700 

0.304403 

9.86982 

8.6403593 

6 

1099 

0. 302277 

9.93970 

8.6485159 

7 

1645 

0.300888 

9.98514 

8.6536345 

8 

2366 

0.299943 

10.0161 

8.6570896 

9 

3290 

0.299274 

10.0380 

8.6595274 

10 

4445 

0.298785 

10.0540 

8. 6613051 

11 

5859 

0.298418 

10.0660 

8.6626398 

Table 

8.2.  The 

values  of 

and  = 

fof’  the  refined  mesh  shown 

Fig.8. 

2 

P 

N 

'^l 

pfq] 

^1 

"“<£(0) 

2 

1438 

0.297669 

10.5183 

8.6049893 

3 

2582 

0.296866 

10. 1308 

8.6320936 

4 

4852 

0.296603 

10. 1317 

8.6671875 

5 

8304 

0.296589 

10. 1210 

8.6686951 

6 

13326 

0.296584 

10. 1214 

8. 6690024 
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We  use  the  values  on  (8.1)  to  compute  the  relative  errors 


-*i'  1^5'”  -S' 


In  Figure  8.4  we  show  these  relative  errors  for  the  mesh  shown  in  Figure  2. 

We  see  that  the  errors  of  aind  decrease  with  the  sajne  rate  as  the 

square  of  the  error  in  the  energy  norm  ais  predicted  by  theorems  6.1,  7.1.  We 

-2/3 

obtain  from  estimates  (5.7)  (5.8)  that  the  error  behaves  like  0((N(q))  ). 

This  rate  is  indicated  in  Figure  8.4  too.  We  see  that  the  aisymptotic  analysis 
of  previous  section  is  completely  applicable  in  the  computed  range. 


In  Figure  8.5  we  show  analogous  results  for  the  refined  mesh.  We  see 

that  also  here  the  accxiracy  of  and  increases  as  the  square  of  the 

error  measured  in  the  energy  norm.  We  au'e  in  the  prezisymptotic  range  of  the 

p-version,  which  can  be  interpreted  as  the  h-p  version.  For  large  p,  the 

-2/3 

errors  again  decrease  with  the  rate  0(N(q)  ).  The  corresponding  slope  is 

indicated  in  Fig.  8.5. 
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Figvire  8.5  Relative  error  for  the  mesh  with  geometric  refinement 

We  note  that  we  have  used  here  the  classical  p-version  with  the  same 
degree  p  of  polynomials  in  all  elements.  By  using  low  values  of  p  close 
to  singularities  and  high  values  of  p  away  from  the  singularities  compau'able 
results  czui  be  obtained  with  a  much  smaller  number  of  degrees  of  freedom. 
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